function [surplus,n1, n2] = demandModel3D(m1,m2,sd1,sd2, r12, w1, w2, p1, p2,  N);
PI = 3.142;

% m1 = 5;    mean willingness to pay data price
% m2 = 5;    mean willingness to pay vedio price
% 
% sd1 = 1;  stand deviation of data willings to pay
% sd2 = 1;  stand deviation of vedio willings to pay
% 
% r12 = -0.5;   corelation number between data and vedio
% 
% w1 = 10;    maximum willingnesss to pay for data
% w2 = 10;    maximum willingnesss to pay for vedio

const1 = 2*PI*sd1*sd2*sqrt(1 - r12*r12);
const1 = 1/const1;

const2 = 2*(1 - r12*r12);
const2 = -1/const2;

n1 = 0;
n2 = 0;

surplus = 0;
for x1=1:1:w1
    for x2=1:1:w2
            y1 = x1 - m1;
            y2 = x2 - m2;
            
            part1 = y1*y1/(sd1*sd1);      	
            part2 = 2*r12*y1*y2/(sd1*sd2);
            part3 = y2*y2/(sd2*sd2);
            
            P(x1,x2) = const1*exp(const2*(part1-part2+part3));   
            val = maxVal(x1-p1, x2-p2);
            if( val == 1)
                n1 = n1 + P(x1,x2)*N;
                surplus = surplus + x1-p1;
            elseif (val == 2)
                n2 = n2 + P(x1,x2)*N;
                surplus = surplus + x2-p2;
            end
    end
end

function val = maxVal(v1, v2);
if( v1 > v2 && v1 >= 0)
    val = 1;
elseif( v2 >= v1 && v2 >= 0)
    val = 2;
else
    val = 3;
end
